Advanced Concepts of Clustering in Insurance: Clustering of Marine Losses

1 Primer

This is an example on how to use clustering techniques to analyze insurance data, not a blue print or best practice that can be followed for any other analysis. Any data analysis, in particular in unsupervised learning, relies heavily on domain knowledge and the application at hand. Even for this data set several other reasonable approaches could have been taken given a different set of assumptions and goals.

2 The dataset

This is an CAS datasets. From the official description:

“The univariate dataset was collected by an unknown French private insurer and comprise 1,274 marine losses between the January 2003 and June 2006. The status of the claim (settled or opened) is determined at the end of June 2006.”

THUMBS UP: Who of you has experience in marine insurance business?

# install.packages("CASdatasets", repos = "http://cas.uqam.ca/pub/", type="source")
library(CASdatasets)

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio

seed_var <- 03052021

data(fremarine)

Data preview:

library(data.table)
cat_vars <- c("ShipBrand","ShipHull","Departement")
num_vars <- c("ShipPower","ShipEngNb","ShipEngYear","ShipBuildYear","ShipLength","ShipTonnage","InsuredValue","ClaimPaid")
dat4clustering <- as.data.table(fremarine)
dat4clustering <- dat4clustering[ClaimStatus=="settled",mget(c(num_vars,cat_vars))]
library(DT)
datatable(dat4clustering, filter = 'top', options = list(pageLength = 5, autoWidth = TRUE))

2.1 Story

We put us into the position of a fictional claims analyst that wants to understand the losses better, having the following questions / whishes in mind:

• Being new in the position, I do not have a good overview of recent losses. Previous analyses where focused on single cases rather then the overall portfolio.

• With a segmentation of our recent losses I could the monitor the loss development over time.

• Unfortunately, we are not enough people to review each claim in depth. New claims that do not fit very well to the existing segmentation, e.g., are quite far away from all centroids, could be marked as outliers and trigger a detailed assessment.

• Claims handlers would like to gain efficiency by specialising on certain portfolio segments (large ships vs small etc.)

• My underwriting colleagues would like to understand our losses better but many of them do not have a quantitative background. With a few representative examples (e.g., cluster centroids), I could give them a simplified but complete overview over our recent losses.

2.2 Data set dictionary

The data set is of mixed data type:

sapply(fremarine,FUN=class)

    OccurDate     ReporDate     ShipCateg     ShipBrand     ShipPower     ShipEngNb   ShipEngYear ShipBuildYear      ShipHull    ShipLength 
       "Date"        "Date"      "factor"      "factor"     "integer"     "integer"     "integer"     "integer"      "factor"     "numeric" 
  ShipTonnage  InsuredValue     ClaimPaid   ClaimCharge ClaimRecourse    ClaimCateg   ClaimStatus    Deductible   HeadQuarter   Departement 
    "numeric"     "numeric"     "numeric"     "numeric"     "numeric"      "factor"      "factor"     "numeric"      "factor"      "factor" 

A (somewhat brief) explanation can be found in the manual and is shown below:

Date variables:

• OccurDate: The day of claim occurence
• ReporDate: The day of claim reporting.ShipCategThe category of the insured ship

Factor variables:

• ShipBrand: The brand of the insured ship
• ShipHull: The hull of the insured ship
• Departement: The French region in which the ship is headquartered

Numeric variables:

• ShipPower:The power of the insured ship
• ShipEngNb: The engine number of the insured ship
• ShipEngYear: The engine year of the insured ship
• ShipBuildYear: The building year of the insured ship
• ShipLength: The length of the insured ship
• ShipTonnage: The tonnage of the insured ship
• InsuredValue: The insured value of the insured ship
• ClaimPaid: The paid amount (EUR) of the claim
• ClaimStatus: The status of the claim. Unsettled

2.3 Assumptions

For simplicity, we assume the paid claims includes the final amount paid to the policyholder after deductible, recourse etc. We will also restrict our analysis to settled claims only for which we have complete and final information.

Additionally, we decide not to use the date since:

• Our (fictional) claims colleagues assure that our recent data set is quite homogeneous in time, there should not be strong trends.
• New claims, obviously, have a date outside of the existing range making it hard to segment them otherwise.
• For ClaimCateg there is no documentation, this category will be ignored.

3 Data preparation

For educational reasons, we distinguish two different settings:

1. The data set does only contain numeric variables.
2. The data set contains numeric and categorical variables.

3.1 Numeric features

3.1.1 Pre-processing

CHAT: Which points do we have to consider in our pre-processing?

We remove categorical variables, missings and ensure that we do not introduce any duplicates

dat_num_no_missings = unique(na.omit(dat4clustering[,mget(num_vars)]))

Here we use the Euclidean metric here. Since the numeric variables are measured on different scales (years, counts, EURs, sizes, etc.) we have to center and scale the data:

dat_num_no_missings_scaled = scale(dat_num_no_missings)

3.1.2 Descriptive analysis

First look at the data:

library(ggridges)
library(ggplot2)

data4ridgeplot <- melt(as.data.table(dat_num_no_missings_scaled))
ggplot(data4ridgeplot) +  geom_density_ridges(aes(x = value, y = variable)) + theme_ridges() + theme(legend.position = "none") + xlim(-2.5,2.5) # given a normal distribution ~1% should be outside of this range

Typically, variables with a large variance turn out to be more important differentiators. Also correlations should be checked.

YES OR NO: Do you expect large correlations?

cor_mat <- cor(dat_num_no_missings_scaled, method="spearman")
library(corrplot)
corrplot.mixed(cor_mat,order="hclust",lower.col = "black")

A combination of histograms and correlations can be obtained using a pairs plot (using Pearson correlation):

library(GGally)

Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

dat_num_no_missings_scaled_df <- as.data.frame(dat_num_no_missings_scaled)
ggpairs(dat_num_no_missings_scaled_df[,c("ShipLength","ShipTonnage")], progress=FALSE)

Further pre-processing steps might be very well reasonable lead to different analyses (cf. this article on researchers degrees of freedom for a general account on this topic). For example:

• Ship length and tonnage are (unsurprisingly) highly correlated. One could remove one of the variables or combine both into a new variable, e.g. length per ton (like in life insurance where one often uses BMI instead of height and weight)
• ClaimPaid variable may be divided by Insured Value. Since these variables are not an objective measure of the ship but a result of an insurance process (underwriting and claims assessment), one may also remove them from the data set. In such a setting, one may further analyze if the clusters based on the similarity of the ships reflect different risk groups, which can be a preliminary step to making a supervised model aiming to predict paid claims based on the ship’s covariates.

In any case, domain knowledge is essential in data pre-processing.

3.2 Mixed-type Data

We’ll filter out NAs here and focus only on the clustering part. In practice, you could use mice to complete the data, but it may be a good starting point to first consider only available data.

dat_no_missings = unique(na.omit(dat4clustering[,mget(c(num_vars,cat_vars))]))

As before the numerical variables have to be scaled with scale() as above. Then we re-add the categorical features again.

dat_no_missings_scaled <- as.data.table(scale(dat_no_missings[,mget(num_vars)]))
dat_no_missings_scaled[,ShipBrand:=dat_no_missings$ShipBrand]
dat_no_missings_scaled[,ShipHull:=dat_no_missings$ShipHull]
dat_no_missings_scaled[,Departement:=dat_no_missings$Departement]

4 Clustering algorithms

4.1 K-Means

YES OR NO: Ever worked with k-means before?

4.1.1 Numeric data

We will perform a K-means Cluster Analysis using the flexclust package here. Below we run clustering algorithms repeatedly for different numbers of clusters and returns the minimum within cluster distance solution for each. The sum of within cluster distances which is decreasing with k (since between cluster distances are not counted). Ideally we’d look for an elbow / hockey stick.

library(flexclust)
k_try <-2:15
find_k <- stepFlexclust(dat_num_no_missings_scaled,k_try,nrep=5,verbose=FALSE, seed=seed_var, FUN=cclust)
plot(find_k)

Based on the result and our applications we select 7 clusters.

nr_clusters = 7
cclust_res <-  getModel(find_k,nr_clusters-1)
cclust_res

kcca object of family ‘kmeans’ 

call:
stepFlexclust(x = dat_num_no_missings_scaled, k = 7L, nrep = 5, verbose = FALSE, FUN = cclust, seed = seed_var)

cluster sizes:

  1   2   3   4   5   6   7 
221 112  55  63 190 137 143 

4.1.2 Mixed data types

One can use the clustMixType as in the first part of the course which uses the k-prototypes algorithm. Alternatively, one can use partitioning around medoids, where one computes the distance matrix using Gower’s distance a-priori. However, this approach is rather slow (and memory intensive) and works only for smaller data sets. We will show an example of a clustering based on a pre-computed distance matrix later in these notes.

4.2 Hierarchical dbscan

YES OR NO: Do you expect the result of h-dbscan to be similar to the one for k-means?

Here we are using a variant of DBSCAN, hierarchical DBSCAN. It has the benefit that we only have to specify the minPts parameter. Note that not all points are clustered, cluster zero contains all outliers but will not be a homogeneous set of similar points.

We want a cluster to contain at least 7 points, thus we set minPts=7. By playing around with this parameter, we see that a larger number produces very few clusters which would not be helpful, and a lower number leads to a lot of outliers. We observe that most points fall into cluster 2. Such a result would not be quite practical if we need segments for our claims colleagues to specialize on. However, it does provide us with a few “uncommon” examples and 48 outlier, we’ll investigate at a later stage.

library(dbscan)
hdbscan_res <- hdbscan(dat_num_no_missings_scaled, minPts = 7)
hdbscan_res

HDBSCAN clustering for 921 objects.
Parameters: minPts = 7
The clustering contains 5 cluster(s) and 48 noise points.

  0   1   2   3   4   5 
 48   8 813   9   7  36 

Available fields: cluster, minPts, cluster_scores, membership_prob, outlier_scores, hc

HDBSCAN essentially computes the hierarchy of all DBSCAN* clusterings, meaning for all (reasonable) values of eps , the radius of the epsilon neighborhood. Then it uses a stability-based extraction method to find optimal cuts in the hierarchy, thus producing a flat solution. Depending on a different value for eps, we could get more or less clusters using the dbscan() function from the same package.

plot(hdbscan_res, show_flat=T)

A higher value of eps would result in only 4 clusters where cluster 4 and 5 will be merged. A lower value will quickly lead to a fine segmentation of cluster 2 with lots of new clusters and potentially, also outliers.

The result seems to be of less use for this data - most ships fall into cluster 2, the other categories are sparsely populated.

Similarly as with k-means, data with mixed data types can per se be used with DBSCAN if the pairwise distances are computed a-priori, which is normally not feasible for large data sets.

4.3 Model-based clustering

4.3.1 Gaussian Mixture Model for numerical data

CHAT: When did Carl Friedrich Gauß define the Normal/Gaussian distribution? (answer at the end of this section)

Here we are making the assumption that each point, given the true cluster assignment, is follows a Gaussian distribution. This stochastic framework allows us to determine the clusters via Maximum Likelihood and evaluate the fit using statistical measures, such as the BIC (Bayesian Information Criterion), which is defined as nr_parameters * log(nr_observations) minus 2 times the log-likelihood.

We fit configurations with 1 to 10 clusters, hence max_clusters=10. Ensure that the number of points >> number of clusters. Most other parameters are used to control the algorithm, that we do not investigate in detail. In practice, some sensitivity analysis will quickly show if deviating from the defaults makes a difference or not. Ideally, one at least observes convergence for large values of km_iter and em_iter. If not, the algorithm may be unstable and not fit for the data at hand.

library(ClusterR)

opt_gmm = Optimal_Clusters_GMM(dat_num_no_missings_scaled, max_clusters = 10, 
                               criterion = "BIC", dist_mode = "eucl_dist", 
                               seed_mode = "random_subset", seed= seed_var,
                               km_iter = 10, em_iter = 10, var_floor = 1e-10,
                               plot_data = T)

The fit becomes better if we increase the number of clusters from one, but eventually the improvement is lower than the penalization in the BIC formula. We’ll decide to use 7 clusters since this value is quite close to the minimum and we want to have a rather low number of clusters for our applications.

gmm_res = GMM(dat_num_no_missings_scaled,gaussian_comps  =  7, dist_mode = "eucl_dist", 
              seed_mode = "random_subset", km_iter = 10, em_iter = 10, seed = seed_var)   

We define a wrapper to easily obtain predictions from our Gaussian Mixture Model (GMM)

predictGMM <- function(x,dat) {
  centroids = x$centroids
  cov = x$covariance_matrices
  w = x$weights
  return (predict_GMM(dat,centroids,cov,w))
}

The resulting clusters are more evenly distributed

pred_gmm_res <- predictGMM(gmm_res,dat_num_no_missings_scaled)
summary(factor(pred_gmm_res$cluster_labels))

  0   1   2   3   4   5   6 
 65  22  81   6 316 235 196 

Since we deal with probabilistic models, we can compute the probability that each point x belongs to cluster k. Thus one often speaks of “soft clustering”. As long as no decision is taken, each point belongs to each cluster with a certain probability. Naturally, one assigns each point to the most likely cluster. We observe that the probabilities of the most likely cluster are nicely skewed to the right, indicating a clear cluster assignment for most points.

probab_of_majority_class <- apply(pred_gmm_res$cluster_proba,1,max)
hist(probab_of_majority_class)

YES OR NO: The 7 cluster distributions have independent components?

gmm_res$covariance_matrices

          [,1]         [,2]      [,3]      [,4]       [,5]        [,6]        [,7]         [,8]
[1,] 0.7045233 0.0000000001 0.9119202 1.5981679 1.41460231 1.605640946  0.56296429  0.043127237
[2,] 2.3970460 1.1809796578 1.5502328 0.8463268 1.41567610 3.897214112 10.48325213  1.094608677
[3,] 0.1494162 0.0000000001 0.8652210 1.0080487 0.22676706 0.013508312  0.11334347  1.697885553
[4,] 0.3497369 3.2622415977 0.6673540 0.4346947 0.95265254 0.640867547  0.41243678 58.533182013
[5,] 0.1441920 0.0000000001 1.2042783 1.2350323 0.08981686 0.003363012  0.01527854  0.001554085
[6,] 0.2207778 0.0000000001 0.7247845 0.6685620 0.03142149 0.013553045  0.06685488  0.003018756
[7,] 0.6988894 0.0000000001 0.9686870 0.4843733 0.56035609 0.904982843  0.66450006  0.014623842

What appears like a non-diagonal covariance matrix, is in fact something else: Each cluster distribution lives in an dim(dat_num_no_missings_scaled)[2] dimensional space, since one needs one dimension per variable, and has a diagonal covariance matrix stated in the rows. So the answer is yes, they do have independent components.

Note that the 2nd component of many clusters has a zero variance, thus there is no variation. The 2nd component refers to the 2nd feature, the number of engines.

Gauß defined the Normal distribution in 1809 in his work “Theoria motus corporum coelestium in sectionibus conicis solem ambientium”, which you can surprisingly even buy on Amazon.

4.3.2 Mixed data types

For data sets with numerical and categorical features, we assume a multivariate distribution that is Gaussian for numerical components and follows a Dirichlet Distribution (that is a multivariate extension of the beta distribution) for categorical features.

library(MixAll)

In contrast to the previous setting with only numerical variables, the optimal number of clusters is 5.

ldata = list(dat_no_missings_scaled[,mget(cat_vars)],
             dat_no_missings_scaled[,mget(num_vars)])
lnames = c("categorical_pk_pjk","gaussian_pk_sjk")
set.seed(seed_var)
clustermix_res <- clusterMixedData(ldata, lnames, 
                                   nbCluster = 5:9, 
                                   strategy = clusterFastStrategy(),
                                   criterion = "BIC")
summary(clustermix_res)

**************************************************************
* model name =  categorical_pk_pjk 
* model name =  gaussian_pk_sjk 
* nbSample       =  789 
* nbCluster      =  5 
* lnLikelihood   =  18043.41 
* nbFreeParameter=  349 
* criterion name =  BIC 
* criterion value=  -33758.72 
**************************************************************

There is no need to re-fit the model with this package:

pred_clustermix <- clustermix_res@zi
# pred_clustermix <- clusterPredict(ldata,clustermix_res) # for new data one would need the predict function
table(pred_clustermix)

pred_clustermix
  1   2   3   4   5 
  3  21  37 242 486 

3 minute BREAK. Have a bio break if you need one, get a tea / coffee or simply a glass of water.

5 Cluster evaluation and interpretation

5.1 Cluster evaluation metrics

We’d focus on k-means and GMM here since hdbscan produced essentially only one large cluster. Moreover, cluster 0 cannot be treated as a cluster since it is merely a collection of outliers. To make a somewhat fair comparison one would have to take out all rows belonging to the outlier cluster 0 before calculating the metrics (thereby accepting that these points are really outliers).

5.1.1 Adjusted RandIndex

The RandIndex goes through all pairs of observations and counts which fraction of them are both in the same cluster or both in different ones, in the k-means or the GMM segmentation (even when the clusters are different, for such a pair the segmentation is consistent since both observations are considered similar or un-similar in each segmentation).

Rand Index Vizualization

There is a rather low consistency between the result from GMM and k-Means, as the adjusted RandIndex shows. The different assumptions lead to quite different segmentations.

external_validation(pred_gmm_res$cluster_labels,flexclust::predict(cclust_res),
                    method = "adjusted_rand_index")

[1] 0.2426295

5.1.2 Davies-Bouldin

library(clusterSim)

The definition of the Davies-Bouldin (DB) index of a segmentation is a bit involved. The dispersion (or variation) within cluster \(i\) can be defined as \(S_i = (\sum_{x \in C_i} |x-c_i|^q)^{1/q}\), where \(c_i\) denotes the centroid of the cluster. For \(q=2\) this corresponds to the standard deviation within cluster \(i\) since the centroid is defined as the mean.

A measure of distance between the centroids of two different clusters can be computed by \(R_{i,j} = (\sum_{k=1}^n |c_{i,k}-c_{j,k}|^p)^{1/p}\). For \(p=2\) this is simply the Euclidean distance between the centroids.

The separation of cluster \(i\) from all other cluster can then be computed by \(r_i = max_{j\neq i}\frac{S_i+S_j}{R_{i,j}}\), that is, for the “worst case” pairing with cluster \(i\), the ratio between the sum of within cluster dispersions and the between cluster distance.

The DB index is then defined as the average over all clusters \(DB = \frac{1}{k}\sum_{i=1}^k r_i\).

Consequently, a lower DB implies a better separation: In such a situation, all clusters are rather far apart and tightly centered around their mean.

cclust_db <- index.DB(dat_num_no_missings_scaled, flexclust::predict(cclust_res),
                      p=2, q=2)
cclust_db$DB

[1] 1.63419

There is a higher value for GMM, meaning that k-means provides a better segmentation.

gmm_db <- index.DB(dat_num_no_missings_scaled, pred_gmm_res$cluster_labels, p=2, q=2)
gmm_db$DB

[1] 2.546584

5.2 Vizualization using T-SNE

Visualization of high-dimensional data into 2 (or 3) dimension can be done via linear (Principal Component Analysis, Multidimensional Scaling) or non-linear projection methods. While non-linear methods often achieve better projections, they are computationally more costly and have several hyper parameters.

T-SNE maps a set of points from a high-dimensional space in a lower-dimensional Mapping should preserve the local neighbourhood structure of each point by minimizing Kullback-Leibler divergence between the two distributions

library(Rtsne)

5.2.1 2 dimensions

set.seed(seed_var)
tsne_out <- Rtsne(dat_num_no_missings_scaled,pca=FALSE, perplexity = 30)

Within this (non-linear) 2-dimensional projection, we see that indeed not all clusters from the k-means segmentation are terribly well separated.

data_tsne <- as.data.table(tsne_out$Y)
ggplot(data_tsne,aes(x=V1,y=V2,color=factor(predict(cclust_res)))) + geom_point() + guides(color=guide_legend(title="Cluster")) +
scale_color_brewer(palette = "Set2")

CHAT: Any ideas how to enhance this visual inspection with (quantitative) evaluation metrics?

We can also connect this with intermediate calculation results of the Davies-Bouldin index. Let us consider Cluster 4, for example.The centroid of this cluster is much closer to the centroid of Cluster 7 than to all other cluster centroids.

cclust_db$d[,4] # matrix of distances between centroids or medoids of clusters

       1        2        3        4        5        6        7 
5.851553 5.895618 6.028385 0.000000 5.953955 6.354263 3.819466 

Some points are far away from the bulk (and close to Cluster 3), hence the dispersion of this cluster is rather high compared to the other clusters (actually highest).

cclust_db$S

[1] 1.133643 1.344648 2.834189 4.761197 1.132800 1.450058 1.427963

Looking at the R matrix, Cluster 4 is indeed not so well separated from Cluster 3 and 7.

cclust_db$R

          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
[1,]       Inf 1.1794494 0.9814852 1.0073976 1.5377024 1.8778583 1.1487073
[2,] 1.1794494       Inf 0.9367326 1.0356583 0.9044690 1.6451230 0.9420658
[3,] 0.9814852 0.9367326       Inf 1.2599372 1.0072547 0.9587183 0.9944728
[4,] 1.0073976 1.0356583 1.2599372       Inf 0.9899299 0.9774942 1.6204255
[5,] 1.5377024 0.9044690 1.0072547 0.9899299       Inf 0.9610392 0.9275255
[6,] 1.8778583 1.6451230 0.9587183 0.9774942 0.9610392       Inf 0.9867794
[7,] 1.1487073 0.9420658 0.9944728 1.6204255 0.9275255 0.9867794       Inf

All in all, in the worst case is aggregation (max of R matrix) is has a lower index than Cluster 1 and 6, for example, which visually show a rather high dispersion and are close to other clusters.

cclust_db$r

[1] 1.877858 1.645123 1.259937 1.620425 1.537702 1.877858 1.620425

With descriptive analyses one can (e.g. histograms and correlations) one can better understand the clusters and their similarity. A bad separation of some clusters can also motivate a re-clustering with a lower number of clusters (investigate clusters with an \(r\)-value that is much higher than the average, i.e., the DB index.)

Of course, we can do the same plot for GMM. We only have to define a helper function the maps the cluster object and a data table to a vector with integers, the cluster assignments per row of the data table.

predictGMM_labels <- function(gmm,newdata){
  return (predictGMM(gmm,newdata)$cluster_labels)
}

data_tsne <- as.data.table(tsne_out$Y)
ggplot(data_tsne,aes(x=V1,y=V2,color=factor(predictGMM_labels(gmm_res,dat_num_no_missings_scaled)))) + geom_point() + guides(color=guide_legend(title="Cluster")) +
  scale_color_brewer(palette = "Set3")

5.2.2 3 dimensions

set.seed(seed_var)
tsne3 <- Rtsne(dat_num_no_missings_scaled,pca=FALSE, perplexity = 30, dims = 3)
data_tsne3 <- as.data.table(tsne3$Y)
data_tsne3$cluster <- factor(predict(cclust_res))

Within this (non-linear) 2-dimensional projection, we see that indeed not all clusters from the k-means segmentation are terribly well separated.

library(plotly)

We’ll now make use of the plotly library, since additional interactivity like rotation and filtering is helpful in a 3 dimensional plot

fig <- plot_ly(data_tsne3, x = ~V1, y = ~V2, z = ~V3, color = ~cluster)
fig <- fig %>%  add_trace(type = 'scatter3d', mode='markers', text = ~cluster,hoverinfo = 'text')
fig

`arrange_()` was deprecated in dplyr 0.7.0.
Please use `arrange()` instead.
See vignette('programming') for more help

5.3 Feature Importance

CHAT: What would you expect to be the most relevant variable? Meaning that if we randomly permute this variable, the clustering changes the most.

Permutation based approach to determine the overall relevance of a feature for the segmentation. Parallels the approach taken in classification tasks.

Source:https://www.modulos.ai/permutation-feature-importance-deep-dive/, retrieved on November 23, 2021

library(FeatureImpCluster)

We observe that build and engine year have the highest importance, while claims paid is irrelevant for the clustering result.

set.seed(seed_var)
FeatureImp_cclust <- FeatureImpCluster(cclust_res,as.data.table(dat_num_no_missings_scaled))
plot(FeatureImp_cclust)

It is very interesting to see the the clusters from GMM rely on quite different variables. Here Length and Tonnage are most relevant, also InsuredValue and ClaimPaid define the segments. On the other hand, build and engine year have a very low relevance.

# Note: this code requires FeatureImpCluster version >= 0.1.5
set.seed(seed_var)
FeatureImp_gmm <- FeatureImpCluster(clusterObj = gmm_res,
                                    data = as.data.table(dat_num_no_missings_scaled),
                                    predFUN = predictGMM_labels)
plot(FeatureImp_gmm)

5.4 Interpretation using rules-based ML models

A simple supervised machine learning algorithm may be easier to interpret than a complex unsupervised algorithm. We train this model on the cluster labels.

CHAT: Which supervised ML models do you know that you would call interpretable?

C5.0 can create an initial tree model then decompose the tree structure into a set of mutually exclusive rules. These rules can then be pruned and modified into a smaller set of potentially overlapping rules. The rules can be created using the rules option.

library(C50)

In the context of machine learning, a rule is a conditional logical statement that can be attached to a predicted value. Note that the conditions have to be considered in the printed order.

set.seed(seed_var)
gmm_rule <- C5.0(x = dat_num_no_missings_scaled, 
                 y = factor(pred_gmm_res$cluster_labels), rules = TRUE)
summary(gmm_rule)

Call:
C5.0.default(x = dat_num_no_missings_scaled, y = factor(pred_gmm_res$cluster_labels), rules = TRUE)


C5.0 [Release 2.07 GPL Edition]     Fri Apr 15 14:59:45 2022
-------------------------------

Class specified by attribute `outcome'

Read 921 cases (9 attributes) from undefined.data

Rules:

Rule 1: (65, lift 14.0)
    ShipPower <= 2.174239
    ShipEngNb > -0.2844257
    ClaimPaid <= 1.552014
    ->  class 0  [0.985]

Rule 2: (11, lift 38.6)
    ShipPower > 3.393784
    ->  class 1  [0.923]

Rule 3: (5, lift 35.9)
    ShipEngNb <= -0.2844257
    ShipLength > 2.181134
    ClaimPaid > 0.5804884
    ->  class 1  [0.857]

Rule 4: (3, lift 33.5)
    ShipPower > 0.3122873
    ShipTonnage <= -0.1037781
    InsuredValue > 0.390325
    ->  class 1  [0.800]

Rule 5: (2, lift 31.4)
    ShipPower > 2.174239
    ShipEngNb > -0.2844257
    ->  class 1  [0.750]

Rule 6: (57, lift 11.2)
    ShipPower <= 0.7142217
    ShipEngNb <= -0.2844257
    ShipTonnage <= -0.1037781
    ClaimPaid > 0.07467122
    ->  class 2  [0.983]

Rule 7: (47, lift 11.1)
    ShipPower <= 0.7142217
    ShipEngNb <= -0.2844257
    ShipTonnage <= -0.3988089
    ClaimPaid > -0.008888884
    ->  class 2  [0.980]

Rule 8: (7, lift 10.1)
    ShipLength <= -0.5353749
    InsuredValue > 0.01789716
    ->  class 2  [0.889]

Rule 9: (8/1, lift 9.1)
    ShipLength > 0.3814468
    ShipTonnage <= -0.2762257
    ->  class 2  [0.800]

Rule 10: (4, lift 127.9)
    ShipLength <= 2.181134
    ShipTonnage > -0.1037781
    ClaimPaid > 0.5804884
    ->  class 3  [0.833]

Rule 11: (2, lift 115.1)
    ShipEngNb > -0.2844257
    ClaimPaid > 1.552014
    ->  class 3  [0.750]

Rule 12: (229/1, lift 2.9)
    ShipPower <= -0.1995777
    ShipEngNb <= -0.2844257
    ShipLength <= -0.161855
    InsuredValue <= -0.3455804
    ClaimPaid <= -0.09244899
    ->  class 4  [0.991]

Rule 13: (232/2, lift 2.9)
    ShipPower <= -0.130871
    ShipEngNb <= -0.2844257
    ShipLength <= -0.5353749
    InsuredValue <= 0.01789716
    ClaimPaid <= -0.008888884
    ->  class 4  [0.987]

Rule 14: (281/3, lift 2.9)
    ShipEngNb <= -0.2844257
    ShipLength <= -0.4335058
    InsuredValue <= -0.2226395
    ClaimPaid <= -0.008888884
    ->  class 4  [0.986]

Rule 15: (194/2, lift 2.9)
    ShipEngNb <= -0.2844257
    ShipLength <= -0.161855
    ShipTonnage <= -0.4652948
    InsuredValue <= -0.3455804
    ClaimPaid <= -0.09244899
    ->  class 4  [0.985]

Rule 16: (177/1, lift 3.9)
    ShipEngNb <= -0.2844257
    ShipLength > -0.4335058
    ShipLength <= 0.3814468
    InsuredValue > -0.3455804
    InsuredValue <= 0.390325
    ClaimPaid <= -0.008888884
    ->  class 5  [0.989]

Rule 17: (19, lift 3.7)
    ShipEngNb <= -0.2844257
    ShipLength > -0.161855
    ShipTonnage <= -0.1037781
    InsuredValue <= -0.3455804
    ClaimPaid <= -0.008888884
    ->  class 5  [0.952]

Rule 18: (8/1, lift 3.1)
    ShipEngNb <= -0.2844257
    ShipLength <= 0.04188319
    ShipTonnage > -0.1037781
    InsuredValue <= 0.8661773
    ->  class 5  [0.800]

Rule 19: (659/441, lift 1.3)
    ShipLength <= 0.04188319
    ->  class 5  [0.331]

Rule 20: (201/10, lift 4.4)
    ShipPower <= 3.393784
    ShipEngNb <= -0.2844257
    ShipTonnage > -0.1037781
    ClaimPaid <= 0.5804884
    ->  class 6  [0.946]

Rule 21: (139/32, lift 3.6)
    ShipPower > 0.7142217
    ShipEngNb <= -0.2844257
    ->  class 6  [0.766]

Default class: 6


Evaluation on training data (921 cases):

            Rules     
      ----------------
        No      Errors

        21   14( 1.5%)   <<


       (a)   (b)   (c)   (d)   (e)   (f)   (g)    <-classified as
      ----  ----  ----  ----  ----  ----  ----
        65                                        (a): class 0
              19                             3    (b): class 1
                    78           2     1          (c): class 2
                           6                      (d): class 3
                               314     2          (e): class 4
                     1           4   230          (f): class 5
                                       1   195    (g): class 6


    Attribute usage:

     95.01% ShipEngNb
     93.38% ClaimPaid
     75.14% ShipLength
     69.82% ShipPower
     57.87% InsuredValue
     53.64% ShipTonnage


Time: 0.0 secs

Interestingly, the attribute usage provides another source to measure variable importance. The lift values can be interpreted as an importance of the rule.

The evaluation of our algorithm should be taken with case, since we are using the training set. The performance can be much worse on a table of new ships. On our trainng set, only 14 ships are not fitted to the cluster they were originally assigned to.

Such an analysis is not limited to C.50, one could use decision tree or a GLM as well.

6 Dimensionality reduction

For production or interpretation, reducing the the problem to a lower dimensional space can be helpful. What does this mean in a cluster analysis? We will show two (out of many) options:

1. Feature Clustering: One can cluster the features to obtain groups of similar features
2. Feature Selection: One can select a subset of features that are most important for the clustering assignment

In any case, this allows to re-compute the clusters on a lower dimensional space.

6.1 Feature Clustering

We start by first computing the distance between the features, therefore we need to transpose the data matrix.

dist_dat <- dist(t(dat_num_no_missings_scaled), method = "euclidean")
dist_dat

              ShipPower ShipEngNb ShipEngYear ShipBuildYear ShipLength ShipTonnage InsuredValue
ShipEngNb      43.08047                                                                        
ShipEngYear    42.50162  40.96290                                                              
ShipBuildYear  39.08921  39.16008    36.65984                                                  
ShipLength     20.69804  39.67558    43.01744      40.79808                                    
ShipTonnage    22.27648  40.80277    43.31603      40.70971   15.27943                         
InsuredValue   28.93815  41.50811    42.26273      38.98007   28.44953    28.35563             
ClaimPaid      41.55550  41.73910    42.89806      42.82165   40.15746    41.09929     42.06968

Based on this distance matrix, we apply a hierarchical clustering.

hc_res <- hclust(dist_dat)
plot(hc_res)

Sine we normalized our data, all features are on the same scale, so we could aggregate features by taking the mean, for example. Thus, in order to reduce the dimension from 8 to 7, we could average Length and Tonnage to a single variable that we call ShipSize.

library(dplyr)

dat_num_aggregated <- data.frame(dat_num_no_missings_scaled)
dat_num_aggregated$ShipSize = 
  (dat_num_aggregated$ShipLength + dat_num_aggregated$ShipTonnage)/2
dat_num_aggregated <- dat_num_aggregated %>% select(!c(ShipLength,ShipTonnage))

Note that the values of this variable can not be interpreted anymore in a meaningful way. In principle, this would also allow us to cluster a ship if one of the values is missing, in this case we would simply put full weight on the other variable.

YES OR NO: Do you expect the new segmentation based on Size instead of Length & Tonnage to be quite similar to the old one such that the RandIndex comparing both segmentations is near to one, say > .9 ?

set.seed(seed_var)
cclust_org_res <- cclust(dat_num_no_missings_scaled,k=7) # eliminate impact of the random seed
set.seed(seed_var)
cclust_aggr_res <- cclust(dat_num_aggregated,k=7)

Despite a negligible reduction (improvement) in the DB index.

cclust_aggr_db <- index.DB(dat_num_aggregated, flexclust::predict(cclust_aggr_res),
                      p=2, q=2)
cclust_aggr_db$DB

[1] 1.575758

cclust_org_db <- index.DB(dat_num_no_missings_scaled, flexclust::predict(cclust_org_res),
                      p=2, q=2)
cclust_org_db$DB

[1] 1.57296

However, the segmentation differs quite a bit from the previous one based on the raw data.

external_validation(flexclust::predict(cclust_org_res),
                    flexclust::predict(cclust_aggr_res),
                    method = "adjusted_rand_index")

[1] 0.7624793

6.2 Feature Selection using Feature Importance

We used Feature Importance as a tool to better understand our clustering results. Here we benefit from it in another way: To un-select features which to not materially contribute to the clustering result.

Iterating the Feature Importance algorithm over various seeds and bootstrap distributions gives a more robust estimate of feature importance.

# randomly sample starting seeds
set.seed(seed_var)
nr_seeds <- 5
sub <- 0.7 # 70% subset in each iteration
biter <- 5 # 5 bootstrap iterations
seeds_vec <- sample(1:10000,nr_seeds)

savedImp <- data.frame(matrix(0,nr_seeds,dim(dat_num_no_missings_scaled)[2]))
count <- 1
for (s in seeds_vec) {
  set.seed(s)
  res <- cclust(dat_num_no_missings_scaled,k=nr_clusters)
  set.seed(s)
  FeatureImp_res <- FeatureImpCluster(res,as.data.table(dat_num_no_missings_scaled),sub = sub,biter = biter)
  savedImp[count,] <- FeatureImp_res$featureImp[sort(names(FeatureImp_res$featureImp))]
  count <- count + 1
}
names(savedImp) <- sort(names(FeatureImp_res$featureImp))

Turns out that feature importance is quite stable. Claims paid might be removed from the data for clustering purposes.

boxplot(savedImp)

set.seed(seed_var)
cclust_noClaimPaid_res <- cclust(data.frame(dat_num_no_missings_scaled) %>% 
                                   select(-ClaimPaid),k=7)

Again there is some difference in the resulting segmentation. It would make sense to further investiage which rows (ships) are assigned to different clusters.

external_validation(flexclust::predict(cclust_org_res),
                    flexclust::predict(cclust_noClaimPaid_res),
                    method = "adjusted_rand_index")

[1] 0.733007

6.3 Feature Selection in Model-based clustering

When fitting our Gaussian Mixture Model (GMM), we used BIC as a measure to decide which number of clusters to use (remember: higher is better). We can do the same here for variable selection.

We’ll use the clustvarsel which offers a greedy search in two directions, starting from a clustering based on a single variable (forward) or from a clustering using all variables (backward), which is what we have done before.

Greedy variable search - general algorithm:

1. Forward initialization: Build the best univariate cluster model: Fit all GMMs with only a single variable, keep the variable \(x_{f_1}\) with the highest BIC.
2. Backward initialization: Fit a GMM based on all variables \((x_{1},\ldots,x_{p})\)
3. Given a fitted GMM computed on the \(m\) variables \((x_{f_1},\ldots,x_{f_m})\), perform an add and remove step and keep the model with the highest BIC, i.e., repeat the following steps:
   • Add step to determine the best larger model: Fit a GMM for all pairs \((x_{f_1},\ldots,x_{f_m},x_j)\) for all \(j\notin \{f_1,\ldots,f_m\}\) and keep the variable \(f_{m+1}:=j\) so that the new model GMM_Add has the highest BIC among all alternatives that include variables \(\{f_1,\ldots,f_m\}\).
   • Remove step to determine the best smaller model: Fit a GMM for all subsets of \(\{f_1,\ldots,f_m\}\) of size \(m-1\), and keep the model GMM_Remove with the highest BIC.
   • Decision: Among GMM, GMM_Add and GMM_Remove, return the model GMM_new with the highest BIC.
   • Stopping criterion: If GMM_new==GMM, stop the search since no improvement is possible.

clustvarsel allows several choices regarding the covariance structure of the Gaussian model, we simply use a covariance matrix that is proportional to the identity matrix here, i.e., an equal variance for each independent component. Further details can be found in the original paper. Of course, one can also compare BIC among different choices to improve model fit.

library(clustvarsel)
library(doParallel)
set.seed(seed_var)
clustvarsel_res_forward <- clustvarsel(dat_num_no_missings_scaled,G=nr_clusters,
                               search="greedy", direction=c("forward"),
                               emModels1 = "E", emModels2 = "EII", # ?mclustModelNames
                               samp = TRUE, parallel = 3)
clustvarsel_res_backward <- clustvarsel(dat_num_no_missings_scaled,G=nr_clusters,
                               search="greedy", direction=c("backward"),
                               emModels1 = "E", emModels2 = "EII", # ?mclustModelNames
                               samp = TRUE, parallel = 3)
save(clustvarsel_res_forward,clustvarsel_res_backward,file="Clustering_2\\FeatureSelectionGMM.Rdata")

Due to the rather long computation time we saved the results a-priori and load them from a file - feel free to run the computation yourself.

load("FeatureSelectionGMM.Rdata")

We can see the trajectory of the search: After having selected two variale, neither add nor remove improve the model.

clustvarsel_res_forward$steps.info

Thus the final clustering is only based on two variables. Note that these are the two most relevant variables in the C5.0 model approximating the cluster result.

clustvarsel_res_forward$subset

ClaimPaid ShipEngNb 
        8         2 

The procedure removes four variables before it stops. Note that step 4 seems redundant. This is because the implementation seems slightly different than explained above: The algorithm iterates between Add and Remove, unless one of the steps yields a model with only a single variable or the full model (thus two adds/removes in the beginning before the iteration starts). It stops when both steps where rejected.

clustvarsel_res_backward$steps.info

Thus the model ends up with four variables, extending the forward model with insured vaulue and ship power.

clustvarsel_res_backward$subset

   ShipPower    ShipEngNb InsuredValue    ClaimPaid 
           1            2            7            8 

Note that this is not well in line with the permutation feature approach, where length and tonnage seemed to be the most relevant variables. Since feature importance depends on the model and the definition of itself, results may often be quite contradicting.

Further side note: Despite the correction with the number of parameters in the BIC formula, overfitting might be an issue. Essentially, each step is a statistical test and it is well know that a large number of joint tests increase the false discovery rate if thresholds (e.g. p-values) are not adjusted.

7 Outlier detection

CHAT: In a single sentence, how would you define the term “outlier”?

7.1 Distance based (k-means)

We’ll compute the distance of each data point to its cluster centroid to identify the top 5 points that are farthest away from it, and in that sense “special”.

centers_matrix <- cclust_res@centers[predict(cclust_res,dat_num_no_missings_scaled), ]
distances <- sqrt(rowSums((dat_num_no_missings_scaled - centers_matrix)^2))
outliers <- order(distances, decreasing=T)[1:5]

print(outliers)

[1] 253 528 634 215  76

For example, data point 253 lies in cluster 4 but has a rather low insured value and a very high paid claim. That is something our claims analyst might want to analyze further.

dat_num_no_missings_scaled[outliers,]

      ShipPower  ShipEngNb ShipEngYear ShipBuildYear  ShipLength ShipTonnage InsuredValue  ClaimPaid
[1,]  1.0611906 -0.2844257  -1.0604693    0.02102725  1.80761377   1.7931457   -0.6761163 21.1685837
[2,] -0.5018873  3.5120392   0.8916138    1.02115421  2.75839177   0.3844772    1.0865760 14.6308412
[3,]  0.8859885 -0.2844257  -1.0604693   -0.16081402 -0.05998588  -0.4154304   10.3404623  1.1966049
[4,]  0.8859885 -0.2844257  -1.0604693   -0.16081402 -0.05998588  -0.4154304   10.3404623 -0.2016342
[5,]  0.8859885 -0.2844257  -1.0604693   -0.16081402 -0.05998588  -0.4154304   10.3404623 -0.1860363

predict(cclust_res,dat_num_no_missings_scaled[outliers,])

[1] 4 3 4 4 4

In our T-SNE map we can highligh in which areas our outlier are lying.

data_tsne$distances <- distances
ggplot(data_tsne,aes(x=V1,y=V2,color=distances)) + geom_point() + scale_color_gradient2()

Note that various other analysis can be made. For example, the share of outliers by cluster.

Also, the threshold “top 5” was completely arbitrary. If there are known outliers, this threshold can be calibrated in a way similar to how this is done for classification models (e.g. considering false positives / negatives on a validation set).

7.2 Using DBSCAN

CHAT: How many calculation are required to obtain the outliers in the DBSCAN segmentation ? (positive number)

DBSCAN is the only algorithm among the ones we are focusing here that has a build-in concept of outlier points. Hence no calculations are required.

outlier_dbscan <- which(hdbscan_res$cluster==0)
outlier_dbscan

 [1]  33  35  47  49  54  58  68  76  83  95 100 116 117 128 173 215 233 240 253 263 300 311 312 321 391 397 422 436 439 505 517 528 586 612
[35] 620 623 634 637 661 750 757 771 774 780 814 863 888 893

All 5 outliers from k-means are also outliers according to the DBSCAN methodology.

outliers %in% outlier_dbscan

[1] TRUE TRUE TRUE TRUE TRUE

7.3 Probabily-based (Model-based clustering)

Definition of outlier: probability of belonging to the assigned cluster is below 50% (this can happen since we have more than 2 clusters)

outlier_gmm <- which(probab_of_majority_class<.5)
outlier_gmm

[1] 321 479 655 658 666 674 705 730 859

The outliers in GMM do not include the 5 outliers from k-means

outliers %in% outlier_gmm

[1] FALSE FALSE FALSE FALSE FALSE

And only the first one is also an outlier in DBSCAN.

outlier_gmm %in% outlier_dbscan

[1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

All in all, we see that there is no universal concept of an outlier. It strongly depends on the clsutering methodology and how outlier is defined based on that.

Ideally, and unsupervised outlier detection approach is used to identify observations for manual inspection. The result of this inspection is tracked so that there will be, in the future, a data base to train a supervised learning algorithm.